---
title: "Sex & statistical power - simulations"
author: "Szymek Drobniak"
date: "`r Sys.Date()`"
format:
html:
toc: true
toc-location: left
toc-depth: 2
theme: simplex
embed-resources: true
code-fold: show
code-tools: true
number-sections: true
crossref:
fig-title: Figure # (default is "Figure")
tbl-title: Table # (default is "Table")
title-delim: — # (default is ":")
fig-prefix: Fig. # (default is "Figure")
tbl-prefix: Tab. # (default is "Table")
editor_options:
chunk_output_type: console
evaluate:
cache: true
---
```{r}
library (ggplot2)
library (pander)
library (emmeans)
library (tidyverse)
library (here)
library (ggpubr)
library (ggridges)
```
Analysis preps (colours, options, auxiliary functions)
```{r}
cbb_palette <- c (
"#000000" , "#E69F00" , "#56B4E9" , "#009E73" , "#F0E442" ,
"#0072B2" , "#D55E00" , "#CC79A7"
)
## Avoid table wrapping
options (width = 800 )
set.seed (999 )
source (here ("PhillipsKarpEtal_data/funs.r" ))
```
Below results repeat simulations from the paper of *Phillips et al. (2018) PLoS Biol 21(6): e3002129*. Original results simplified by removing the post-hoc and the (pooled) t-test methods of analysing data. Generation of simulated data was modified to use a formal linear model defined as:
$$\mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e}$$
with $\mathbf{b}=(\mu,\beta_{Treat},\beta_{Sex},\beta_{Treat:Sex})$ and $\mathbf{e} \sim \mathcal{N}(\mathbf{0}, \mathbf{s})$, with $\mathbf{s}=\sigma\sqrt{\mathbf{1}+\mathbf{X}_{Sex}\cdot\alpha}$ ($\alpha$ quantifies the amount of heteroscedasticity in sex variances, with **SD** $\sigma_{m}=\sigma_{f}\sqrt{1+\alpha}$, or male **variance** being $\alpha\cdot100\%$ greater than female). Regression coefficients ($\mathbf{b}$) are compatible with contrasts (sex/treatment effects) used in Phillips et al. (2023).
# Scenarios from Fig. 3
## From original paper
**Assuming homoscedasticity of sex variances**
::: {.panel-tabset}
## No sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + 0
all_effs_sex_none <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none <- all_effs_sex_none %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("none" ))
calc_power_sex_none %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1
all_effs_sex_small <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small <- all_effs_sex_small %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("none" ))
calc_power_sex_small %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Large sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1
all_effs_sex_large <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large <- all_effs_sex_large %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("none" ))
calc_power_sex_large %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Plots combining all scenarios from Fig. 1
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_none, calc_power_sex_small,
calc_power_sex_large)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (sex_eff_size = factor (sex_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = sex_eff_size,
group = interaction (sex_eff_size, method))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = method)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Sex effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Stat method" ) + ylim (0 , 1 )
all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate (sex_eff_size = factor (sex_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect,
group = interaction (sex_eff_size, Effect))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = sex_eff_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Sex effect size" ) + ylim (0 , 1 )
```
## Heteroscedastic sex variances
We assume several scenarios with $\alpha$ equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
**Low heteroscedasticity**
```{r}
alpha <- 1
```
::: {.panel-tabset}
## No sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_none_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none_lowh <- all_effs_sex_none_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_none_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_small_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small_lowh <- all_effs_sex_small_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_small_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Large sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_large_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
**Large heteroscedasticity**
```{r}
alpha <- 2
```
::: {.panel-tabset}
## No sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_none_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none_largeh <- all_effs_sex_none_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_none_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_small_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small_largeh <- all_effs_sex_small_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_small_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Large sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_large_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
It should be of note that evolutionary biology has seen heteroscedasticities exceeding 100% in variance ration between sexes. Combined results below show the impact of heteroscedasticity.
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_small,
calc_power_sex_large,
calc_power_sex_small_lowh,
calc_power_sex_large_lowh,
calc_power_sex_small_largeh,
calc_power_sex_large_largeh)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (sex_eff_size = factor (sex_eff_size,
levels = c ("small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = sex_eff_size,
group = interaction (sex_eff_size, heterosc))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Sex effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
plotA <- all_sex_effs_sex_only %>%
filter (sex_eff_size %in% c ("large" )) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
plotA
```
# Scenarios from Fig. 4 - interaction of magnitude
## From original paper
Below simulations were only slightly modified. We assumed a moderate (+0.5) overall effect of sex and a varying interaction increasing the magnitude of sex difference in the treated group (+0, +0.5, +1).
**Assuming homoscedasticity of sex variances**
::: {.panel-tabset}
## No sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + 0
all_effs_sex_ix_none <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1
all_effs_sex_ix_small <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Large sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1
all_effs_sex_ix_large <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Plots combining all scenarios from Fig. 1
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_ix_none, calc_power_sex_ix_small,
calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction (ix_eff_size, method))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = method)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Ix effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Stat method" ) + ylim (0 , 1 )
all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect,
group = interaction (ix_eff_size, Effect))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = ix_eff_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Ix effect size" ) + ylim (0 , 1 )
```
## Heteroscedastic sex variances
We assume several scenarios with $\alpha$ equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
**Low heteroscedasticity**
```{r}
alpha <- 1
```
::: {.panel-tabset}
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_ix_none_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, small ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1 + alpha
all_effs_sex_ix_small_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, large ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1 + alpha
all_effs_sex_ix_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
**Large heteroscedasticity**
```{r}
alpha <- 2
```
::: {.panel-tabset}
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_ix_none_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, small ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1 + alpha
all_effs_sex_ix_small_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, large ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1 + alpha
all_effs_sex_ix_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Power consequences of heteroscedasticities with varying (magnitude-only) interaction.
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_ix_small,
calc_power_sex_ix_large,
calc_power_sex_ix_small_lowh,
calc_power_sex_ix_large_lowh,
calc_power_sex_ix_small_largeh,
calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction (ix_eff_size, heterosc))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Ix effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
plotB <- all_sex_effs_sex_only %>%
filter (ix_eff_size %in% c ("large" )) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
plotB
```
# Scenarios from Fig. 5 - interaction of direction
## From original paper
Below simulations were only slightly modified. We assumed a strong (-1) overall effect of sex and a varying interaction resulting of crossing of sex effects on the treatment effect (0, -1, -2).
**Assuming homoscedasticity of sex variances**
::: {.panel-tabset}
## No sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + 0
all_effs_sex_ix_none <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 1
sex_sd <- 1
all_effs_sex_ix_small <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Large sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 2
sex_sd <- 1
all_effs_sex_ix_large <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("none" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Plots combining all scenarios from Fig. 1
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_ix_none, calc_power_sex_ix_small,
calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction (ix_eff_size, method))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = method)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Ix effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Stat method" ) + ylim (0 , 1 )
all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("none" , "small" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect,
group = interaction (ix_eff_size, Effect))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = ix_eff_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Ix effect size" ) + ylim (0 , 1 )
```
## Heteroscedastic sex variances
We assume several scenarios with $\alpha$ equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
**Low heteroscedasticity**
```{r}
alpha <- 1
```
::: {.panel-tabset}
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_ix_none_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, small ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 1
sex_sd <- 1 + alpha
all_effs_sex_ix_small_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, large ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 2
sex_sd <- 1 + alpha
all_effs_sex_ix_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("low" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large_lowh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
**Large heteroscedasticity**
```{r}
alpha <- 2
```
::: {.panel-tabset}
## Small sex eff, no ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_ix_none_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" ))
calc_power_sex_ix_none_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, small ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 1
sex_sd <- 1 + alpha
all_effs_sex_ix_small_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("small" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("small" ))
calc_power_sex_ix_small_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## Small sex eff, large ix
```{r}
#| code-fold: true
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 2
sex_sd <- 1 + alpha
all_effs_sex_ix_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (heterosc = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" ))
calc_power_sex_ix_large_largeh %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Power consequences of heteroscedasticities with varying (direction-only) interaction. Power anomalies indicate pretty low inferential abilities of such models at low sample sizes.
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (calc_power_sex_ix_small,
calc_power_sex_ix_large,
calc_power_sex_ix_small_lowh,
calc_power_sex_ix_large_lowh,
calc_power_sex_ix_small_largeh,
calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter (Effect %in% c ("Treatment" )) %>%
mutate (ix_eff_size = factor (ix_eff_size,
levels = c ("small" , "large" ))) %>%
mutate (method = recode (method,
"ANOVA" = "Factorial" )) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction (ix_eff_size, heterosc))) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Ix effect size" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
all_sex_effs_sex_only %>%
filter (ix_eff_size %in% c ("large" )) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate (heterosc = factor (heterosc,
levels = c ("none" , "low" , "large" ))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = heterosc)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Heteroscedasticity" ) + ylim (0 , 1 )
```
# Impact of sample size on power
## Large sex eff, no ix
Below simulation code generates power estimates for the following sample size scenarios: 5M+5F in each treatment group (original scenario, totql sample 20 individuals), 10M+10F (N = 40), 50M+50F (N = 200).
First - low heteroscedasticity:
::: {.panel-tabset}
## 5M+5F
```{r}
#| code-fold: true
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_lowh_5 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_5 <- all_effs_sex_large_lowh_5 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("5M+5F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_large_lowh_5 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 10M+10F
```{r}
#| code-fold: true
no_per_gp <- 10
all_effs_sex_large_lowh_10 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_10 <- all_effs_sex_large_lowh_10 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("10M+10F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_large_lowh_10 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 50M+50F
```{r}
#| code-fold: true
no_per_gp <- 50
all_effs_sex_large_lowh_50 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_50 <- all_effs_sex_large_lowh_50 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("50M+50F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_large_lowh_50 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Large heteroscedasticity:
::: {.panel-tabset}
## 5M+5F
```{r}
#| code-fold: true
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_largeh_5 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_5 <- all_effs_sex_large_largeh_5 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("5M+5F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_large_largeh_5 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 10M+10F
```{r}
#| code-fold: true
no_per_gp <- 10
all_effs_sex_large_largeh_10 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_10 <- all_effs_sex_large_largeh_10 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("10M+10F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_large_largeh_10 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 50M+50F
```{r}
#| code-fold: true
no_per_gp <- 50
all_effs_sex_large_largeh_50 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_50 <- all_effs_sex_large_largeh_50 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (sample_size = rep ("50M+50F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_large_largeh_50 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
## Large sex effect, crossed reaction norms (interaction)
Identical sample size scenarios.
First - low heteroscedasticity:
::: {.panel-tabset}
## 5M+5F
```{r}
#| code-fold: true
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 2
sex_sd <- 1 + alpha
all_effs_sex_ix_large_lowh_5 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_5 <- all_effs_sex_ix_large_lowh_5 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("5M+5F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_ix_large_lowh_5 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 10M+10F
```{r}
no_per_gp <- 10
all_effs_sex_ix_large_lowh_10 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_10 <- all_effs_sex_ix_large_lowh_10 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("10M+10F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_ix_large_lowh_10 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 50M+50F
```{r}
no_per_gp <- 50
all_effs_sex_ix_large_lowh_50 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_50 <- all_effs_sex_ix_large_lowh_50 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("50M+50F" )) %>%
mutate (heterosc = rep ("low" ))
calc_power_sex_ix_large_lowh_50 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
Large heteroscedasticity:
::: {.panel-tabset}
## 5M+5F
```{r}
#| code-fold: true
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- - 2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 2
sex_sd <- 1 + alpha
all_effs_sex_ix_large_largeh_5 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_5 <- all_effs_sex_ix_large_largeh_5 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("5M+5F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_ix_large_largeh_5 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 10M+10F
```{r}
no_per_gp <- 10
all_effs_sex_ix_large_largeh_10 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_10 <- all_effs_sex_ix_large_largeh_10 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("10M+10F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_ix_large_largeh_10 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
## 50M+50F
```{r}
no_per_gp <- 50
all_effs_sex_ix_large_largeh_50 <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_50 <- all_effs_sex_ix_large_largeh_50 %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("large" )) %>%
mutate (sample_size = rep ("50M+50F" )) %>%
mutate (heterosc = rep ("large" ))
calc_power_sex_ix_large_largeh_50 %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (group = Effect)) +
theme_classic () +
ylab ("Power" ) + ylim (0 , 1 )
```
:::
## Summaries for sample size scenarios
No interaction, low (left) and large (right) heteroscedasticity:
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (
calc_power_sex_large_lowh_5,
calc_power_sex_large_lowh_10,
calc_power_sex_large_lowh_50,
calc_power_sex_large_largeh_5,
calc_power_sex_large_largeh_10,
calc_power_sex_large_largeh_50
)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
filter (heterosc %in% c ("low" )) %>%
mutate (sample_size = factor (sample_size,
levels = c ("5M+5F" , "10M+10F" , "50M+50F" ))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = sample_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Sample size" ) + ylim (0 , 1 )
plot2 <- all_sex_effs_sex_only %>%
filter (heterosc %in% c ("large" )) %>%
mutate (sample_size = factor (sample_size,
levels = c ("5M+5F" , "10M+10F" , "50M+50F" ))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = sample_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Sample size" ) + ylim (0 , 1 )
ggarrange (plot1, plot2, ncol = 2 , nrow = 1 , common.legend = TRUE , legend = "right" )
```
Interaction of effect direction present, low (left) and high (right) heteroscedasticity:
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (
calc_power_sex_ix_large_lowh_5,
calc_power_sex_ix_large_lowh_10,
calc_power_sex_ix_large_lowh_50,
calc_power_sex_ix_large_largeh_5,
calc_power_sex_ix_large_largeh_10,
calc_power_sex_ix_large_largeh_50
)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
filter (heterosc %in% c ("low" )) %>%
mutate (sample_size = factor (sample_size,
levels = c ("5M+5F" , "10M+10F" , "50M+50F" ))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = sample_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Sample size" ) + ylim (0 , 1 )
plot2 <- all_sex_effs_sex_only %>%
filter (heterosc %in% c ("large" )) %>%
mutate (sample_size = factor (sample_size,
levels = c ("5M+5F" , "10M+10F" , "50M+50F" ))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = p_value, colour = Effect)) +
geom_hline (yintercept = 0.8 , linetype = "dashed" , colour = "gray" ) +
geom_line (aes (linetype = sample_size)) +
theme_classic () +
ylab ("Statistical power" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) +
scale_linetype_discrete (name = "Sample size" ) + ylim (0 , 1 )
ggarrange (plot1, plot2, ncol = 2 , nrow = 1 , common.legend = TRUE , legend = "right" )
```
# Simulating sample size
## Simulations assuming target power = 80%
Scenario 1: strong sex effect, none to large heteroscedasticity. Assumed power 80% for treatment effect.
```{r}
#| code-fold: true
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_noh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("none" ))
```
```{r}
#| code-fold: true
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("low" ))
```
```{r}
#| code-fold: true
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1 + alpha
all_effs_sex_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("large" ))
```
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Heteroscedasticity" , values = cbb_palette) + coord_trans (y = "log10" ) +
theme (text = element_text (size= 16 ))
plot1
plot1.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Heteroscedasticity" , values = cbb_palette) +
xlim (0.3 , 0.7 ) + ylim (5 , 70 )
plot1.1
```
```{r}
#| code-fold: true
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_noh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("none" ))
```
```{r}
#| code-fold: true
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("low" ))
```
```{r}
#| code-fold: true
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("large" ))
```
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot2 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) + coord_trans (y = "log10" ) +
theme (text = element_text (size= 16 ))
plot2
plot2.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Heteroscedasticity" , values = cbb_palette) +
xlim (0.3 , 0.7 ) + ylim (5 , 25 )
plot2.1
```
```{r}
#| code-fold: true
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_noh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("none" ))
```
```{r}
#| code-fold: true
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_lowh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("low" ))
```
```{r}
#| code-fold: true
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- - 0.5
sex_sd <- 1 + alpha
all_effs_sex_large_largeh <- rbind (
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed (
n_rep = 1000 , no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9 , sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9 , sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter (! Effect == "Residuals " ) %>%
group_by (Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all (.funs = list (~ sum (. < 0.05 , na.rm = TRUE ) / n ())) %>%
mutate (sex_eff_size = rep ("large" )) %>%
mutate (ix_eff_size = rep ("none" )) %>%
mutate (heterosc = rep ("large" ))
```
```{r}
#| code-fold: true
all_sex_effs_sex_only <- rbind (
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame (lapply (all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot3 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Effect" , values = cbb_palette) + coord_trans (y = "log10" ) +
theme (text = element_text (size= 16 ))
plot3
plot3.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot (aes (x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line () +
theme_bw () +
ylab ("Sample size" ) +
xlab ("Main treatment effect" ) +
scale_colour_manual (name = "Heteroscedasticity" , values = cbb_palette) +
xlim (0.3 , 0.7 ) + ylim (5 , 2200 )
plot3.1
```
## Combined plots
(In all below examples the y axis (sample size) indicates the number of individuals per group and sex - hence the total sample size is 4 times the number indicated on the y axis.)
Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for full range of effect sizes and on log-scale for sample size:
```{r}
#| fig-height: 13
comb1 <- ggarrange (plot1, plot2, plot3, ncol = 1 , nrow = 3 , common.legend = TRUE , legend = "bottom" , labels = c ("A" , "B" , "C" ))
comb1
```
Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for moderate-sized effect sizes using identity scale for sample size:
```{r}
#| fig-height: 13
comb2 <- ggarrange (plot1.1 , plot2.1 , plot3.1 , ncol = 1 , nrow = 3 , common.legend = TRUE , legend = "bottom" , labels = c ("A" , "B" , "C" ))
comb2
```
Combined plots demonstrating the impact of heteroscedasticity on power for a no-interaction case (A) and interaction-of-magnitude case (B; stronger effect in the more variable sex):
```{r}
comb3 <- ggarrange (plotA, plotB, ncol = 2 , nrow = 1 , common.legend = TRUE , legend = "bottom" , labels = c ("A" , "B" ))
comb3
```
# Source of heterogeneity
Below plots show hypothetical scenarios when heterogeneity is a natural outcome of the underlying statistical process. Diagonal lines connect averages belonging to one treatment group (control vs. experimental). Thick vertical lines show SD of each subgroup.
## Poisson distribution
Simulation draws from $Pois(\lambda_i)$ assumes sex/treatment specific averages ($\lambda_i$) - see `lambda1...lambda4` below.
```{r}
N <- 500
lambda1 <- 3
lambda2 <- 9
lambda3 <- 8
lambda4 <- 15
toy_poisson <- data.frame (
Treatment = c (rep ("Ctrl" , 2 * N), rep ("Exp" , 2 * N)),
Sex = c (rep ("F" , N), rep ("M" , N), rep ("F" , N), rep ("M" , N)),
Response = c (rpois (N, lambda1), rpois (N, lambda2), rpois (N, lambda3), rpois (N, lambda4))
)
toy_poisson$ Group <- interaction (toy_poisson$ Sex, toy_poisson$ Treatment)
```
```{r}
plot1 <- ggplot (data = toy_poisson, aes (x = Response, y = Group)) +
geom_density_ridges (aes (fill = Sex, scale = 0.8 ),
stat = "binline" , bins = 15 , alpha = 0.6 , draw_baseline = FALSE
) +
theme_ridges () +
geom_point (size = 4 , stat = "summary" , fun = "mean" ) +
stat_summary (
fun = "mean" , geom = "line" , linewidth = 1 ,
aes (group = Treatment, linetype = Treatment)
) +
stat_summary (
fun.max = function (x) mean (x) + sd (x),
fun.min = function (x) mean (x) - sd (x), geom = "errorbar" ,
linewidth = 1.5 , width = 0
) +
theme_bw () +
coord_flip () +
theme (text = element_text (size = 16 ))
plot1
```
## Log-normal distribution
Simulation samples from $Lognormal(\mu_i, \sigma^2=0.1)$ with sex/treatment specific means (`mu1...mu4` ) below.
```{r}
N <- 500
mu1 <- 0.01
mu2 <- 0.5
mu3 <- 0.4
mu4 <- 0.8
sd <- 0.31
toy_lognorm <- data.frame (
Treatment = c (rep ("Ctrl" , 2 * N), rep ("Exp" , 2 * N)),
Sex = c (rep ("F" , N), rep ("M" , N), rep ("F" , N), rep ("M" , N)),
Response = c (rlnorm (N, mu1, sd), rlnorm (N, mu2, sd), rlnorm (N, mu3, sd), rlnorm (N, mu4, sd))
)
toy_lognorm$ Group <- interaction (toy_poisson$ Sex, toy_poisson$ Treatment)
```
```{r}
plot2 <- ggplot (data = toy_lognorm, aes (x = Response, y = Group)) +
geom_density_ridges (aes (fill = Sex, scale = 0.8 ),
stat = "binline" , bins = 15 , alpha = 0.6 , draw_baseline = FALSE
) +
theme_ridges () +
geom_point (size = 4 , stat = "summary" , fun = "mean" ) +
stat_summary (
fun = "mean" , geom = "line" , linewidth = 1 ,
aes (group = Treatment, linetype = Treatment)
) +
stat_summary (
fun.max = function (x) mean (x) + sd (x),
fun.min = function (x) mean (x) - sd (x), geom = "errorbar" ,
linewidth = 1.5 , width = 0
) +
theme_bw () +
coord_flip () +
theme (text = element_text (size = 16 ))
plot2
```
## Combined plot
```{r}
#| fig-label: fig:source_hetero1
#| fig-cap: Simulated sex-by-treatment data using a Poisson (A) and log-normal (B) distribution. Diagonal lines connect averages belonging to one treatment group (control vs. experimental). Thick vertical lines show SD of each subgroup.
ggarrange (plot1, plot2, ncol = 2 , nrow = 1 ,
common.legend = TRUE , legend = "bottom" , labels = c ("A" , "B" ))
```
## Non-linear response in sexes
Heterogeneity in variances may also arise if the underlying phenotypic mapping function is non linear and has saturating/desaturating character in a range modified by sex dimorphism and/or experimental treatment.
Here we assume a logistic function:
$$f(x)=\frac{L}{e^{-k(x-x_0)}}$$
with parameter values: `L = 3.5` , `k = 1` , `x0 = 1` .
```{r}
mapping_logistic <- function (x, L = 3.5 , x0 = 1 , k = 1 ) L/ (1 + exp (- 1 * k* (x - x0)))
```
Using this function, we can generate data assuming some treatment and sex differences:
```{r}
SD_N <- 0.5
mu_env1 <- 1
mu_env2 <- 2
mu_env3 <- 3
mu_env4 <- 4.5
toy_nonlin <- data.frame (
Treatment = c (rep ("Ctrl" , 2 * N), rep ("Exp" , 2 * N)),
Sex = c (rep ("F" , N), rep ("M" , N), rep ("F" , N), rep ("M" , N)),
Env = c (
rnorm (N, mu_env1, SD_N), rnorm (N, mu_env2, SD_N),
rnorm (N, mu_env3, SD_N), rnorm (N, mu_env4, SD_N)
)
)
toy_nonlin$ Response <- mapping_logistic (toy_nonlin$ Env)
```
```{r}
#| fig-label: source_hetero2
#| fig-cap: Heteroscedasticity between sexes (colours) and treatment groups (solid/dashed lines) resulting from a non-linear (here - logistic) mapping of some manipulated environemntal variable into the measured response. Thin dotted line - underlying mapping function. Points - realisations of mapping for simulated environmental data. Histograms relate to observed response measurements.
ggplot () +
geom_point (
data = toy_nonlin,
aes (y = Env, x = Response, colour = Sex),
pch = 1 , size = 3 , alpha = 0.8
) +
geom_line (
colour = "black" , linetype = 2 ,
aes (y = seq (- 1 , 8 , 0.1 ), x = mapping_logistic (seq (- 1 , 8 , 0.1 )))
) +
geom_density_ridges (
data = toy_nonlin %>%
group_by (Sex, Treatment) %>%
mutate (meanEnv = mean (Env)),
aes (x = Response, group = meanEnv, y = meanEnv, fill = Sex),
stat = "binline" , alpha = 0.8 , draw_baseline = F
) +
geom_point (
data = toy_nonlin %>%
group_by (Sex, Treatment) %>%
summarise (meanEnv = mean (Env), meanResp = mean (Response)),
aes (x = meanResp, y = meanEnv),
size = 4
) +
geom_line (
data = toy_nonlin %>%
group_by (Sex, Treatment) %>%
summarise (meanEnv = mean (Env), meanResp = mean (Response)),
aes (x = meanResp, y = meanEnv, linetype = Treatment),
linewidth = 1
) +
geom_errorbarh (
data = toy_nonlin %>%
group_by (Sex, Treatment) %>%
summarise (meanEnv = mean (Env), meanResp = mean (Response),
minResp = mean (Response) - sd (Response),
maxResp = mean (Response) + sd (Response)),
aes (y = meanEnv, xmin = minResp, xmax = maxResp),
linewidth = 1.5 , height = 0
) +
# stat_summary(data = toy_nonlin %>%
# group_by(Sex, Treatment) %>%
# mutate(meanEnv = mean(Env)),
# fun = "mean", geom = "line", linewidth = 1,
# aes(x = Response, group = meanEnv, y = meanEnv, linetype = Treatment)
# ) +
coord_flip () + ylab ("Manipulated environmental variable" ) +
theme_ridges () +
theme_bw () + theme (text = element_text (size = 16 ))
```